#include <stdio.h>
#include <assert.h>
#include <math.h>
#define N 26

double f1(double x,double a,double b)
{
	assert(x >= a && x <= b);
	return 1/x;
}

double f2(double x,double a,double b)
{
	assert(x >= a && x <= b);
	return 4/(1+x*x);
}

double Romberg(double a,double b,int n,double (*f)(double ,double ,double ))
{
	double h = b - a;
	double t[n][n];
	int k,i;
	t[0][0] = h*(f(a,a,b)+f(b,a,b))/2;
	for(k = 1;k < n;k++)
	{
		double q = 0;
		assert(h * pow(2,k-1) == b-a);
		for(i = 0;i < pow(2,k-1);i++)
			q += f(a+h*(i+1/2),a,b);
		t[0][k] = t[0][k-1]/2 + h*q/2;
		for(i = 1;i <= k;i++)
			t[i][k-i] = (pow(4,i)*t[i-1][k-i+1] - t[i-1][k-i])/(pow(4,i)-1);
		h /= 2;
	}
	return t[n-1][0];
}
int main()
{
	double s;
	double t;
	s = Romberg(1,2,N,f1);
	t = log(2.0);
	printf("Romberg ln2：%.16f\n",s);
	printf("Math ln2：%.16f\n",t);
	s = Romberg(0,1,N,f2);
	t = 3.141593;
	printf("Romberg PI：%.16f\n",s);
	printf("Math PI：%.16f\n",t);	
	return 0;
}
